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Abstract —  Speckle  noise  removal  by  means  of  digital  image 
processors  could  improve  the  diagnostic  potential  of  medical 
ultrasound.  This  paper  addresses  the  speckle  suppression 
issue  within  the  framework  of  wavelet  analysis.  As  a  first 
step  of  our  approach,  the  logarithm  of  the  original  image 
is  decomposed  into  several  scales  through  a  multiresolution 
analysis  employing  the  2-D  wavelet  transform.  Then,  we 
design  a  maximum  a  posteriori  (MAP)  estimator,  which  re¬ 
lies  on  a  recently  introduced  statistical  representation  for 
the  wavelet  coefficients  of  ultrasound  images  [1].  We  use  an 
alpha-stable  model  to  develop  a  blind  noise-removal  proces¬ 
sor  that  performs  a  non-linear  operation  on  the  data.  Fi¬ 
nally,  we  compare  our  technique  to  current  state-of-the-art 
denoising  methods  applied  on  actual  ultrasound  images  and 
we  find  it  more  effective,  both  in  terms  of  speckle  reduction 
and  signal  detail  preservation. 

Keywords  -  speckle  noise,  wavelet  transform,  alpha-stable  dis¬ 
tributions,  MAP  estimation 

I.  Introduction 

Nowadays,  clinicians  are  allowed  to  noninvasively  eval¬ 
uate  physiological  processes  within  the  human  body  by 
means  of  various  medical  imaging  modalities,  such  as  com¬ 
puted  tomography,  positron  emission  tomography,  func¬ 
tional  magnetic  resonance  imaging,  and  ultrasonography. 
Among  them,  ultrasound  imaging  is  definitely  the  one  that 
offers  the  best  price-to-performance  ratio  receiving  thus  an 
added  attention.  However,  one  major  issue  when  using 
this  imaging  modality  is  the  inherent  presence  of  speckle 
noise.  Its  occurrence  is  often  undesirable,  since  it  affects 
the  tasks  of  human  interpretation  and  diagnosis.  On  the 
other  hand,  its  texture  carries  important  information  about 
the  tissue  being  imaged.  Speckle  filtering  is  thus  a  criti¬ 
cal  pre-processing  step  in  medical  ultrasound  imagery,  pro¬ 
vided  that  the  features  of  interest  for  diagnosis  are  not  lost. 

Current  speckle  reduction  methods  are  based  on  tempo¬ 
ral  averaging  [2],  [3],  median  filtering  [4],  and  homomor¬ 
phic  Wiener  filtering  [5].  However,  it  is  recognized  that 
standard  noise  filtering  methods  often  result  in  blurred  im¬ 
age  features.  Indeed,  single-scale  representations  of  signals, 
either  in  time  or  in  frequency,  are  often  inadequate  when 
attempting  to  separate  signals  from  noisy  data. 

Recently,  the  wavelet  transform  has  been  proposed  as 
a  useful  processing  tool  for  signal  recovery  [6],  [7],  [8], 
[9].  Current  state-of-the-art  wavelet-based  speckle  reduc¬ 
tion  and  image  enhancement  techniques  employ  a  combi¬ 
nation  of  wavelet  shrinkage  by  soft  and  hard  thresholding 
together  with  a  generalized  adaptive  gain  (GAG)  for  fea¬ 
ture  emphasis  [8] .  Also,  in  [9]  the  authors  describe  a  com- 

A.  Achim  is  sponsored  by  the  State  Scholarship  Foundation  of 
Greece  (IKY)  under  grant  542  /  1999. 

The  work  of  Dr.  Tsakalides  was  supported  by  the  Greek  General 
Secretariat  for  Research  and  Technology  under  Program  EIIET  II, 
Code  97EA  -  152. 


bination  of  adaptive  weighted  median  filtering  [4]  and  soft 
thresholding.  These  methods  try  to  address  the  inability  of 
the  original  soft  thresholding  technique  to  balance  between 
speckle  suppression  and  signal  detail  preservation. 

In  a  recent  work  [1],  we  have  shown  that  a  successful 
ultrasound  imaging  algorithm  can  achieve  both  noise  re¬ 
duction  and  feature  preservation  if  it  takes  into  considera¬ 
tion  the  true  statistics  of  the  signal  and  noise  components. 
Specifically,  we  have  shown  that  the  subband  decomposi¬ 
tions  of  ultrasound  images  have  significantly  non-Gaussian 
statistics  that  are  best  described  by  families  of  heavy-tailed 
distributions  such  as  the  alpha-stable  and  consequently  we 
designed  a  Bayesian  estimator  that  exploits  these  statistics. 

The  approach  presented  here  is  similar  to  the  method  re¬ 
ported  in  [1].  The  differences  are  that:  (i)  we  select  a  differ¬ 
ent  cost  function  for  the  design  of  the  Bayesian  processor, 
which  gives  raise  to  slightly  different  shapes  of  the  nonlin¬ 
earities  applied  to  the  noisy  wavelet  coefficients,  and  (ii)  we 
propose  a  new  method  for  estimating  the  parameters  of  the 
alpha-stable  distribution  from  noisy  observations,  which  is 
based  on  Koutrouvelis’  [10]  regression  method. 

II.  Preliminaries 

Parametric  Bayesian  processing  presupposes  proper 
modeling  for  the  prior  probability  density  function  (PDF) 
of  the  signal  and  noise.  The  statistical  properties  of  speckle 
noise  were  studied  by  Goodman  [2].  He  has  shown  that,  if 
the  number  of  scatterers  per  resolution  cell  is  large,  a  fully 
developed  speckle  pattern  can  be  modeled  as  the  magni¬ 
tude  of  a  complex  Gaussian  held  with  independent  and 
identically  distributed  (i.i.cl.)  real  and  imaginary  compo¬ 
nents.  In  general,  speckle  noise  has  a  spatial  correlation 
length  on  each  axis  which  is  roughly  the  same  as  the  reso¬ 
lution  cell  size  [11].  In  order  to  generate  spatially  correlated 
speckle  noise  for  use  in  simulations,  one  can  lowpass  filter 
a  complex  Gaussian  random  held  and  take  the  magnitude 
of  the  hltered  output  [12],  [13]. 

Arsenault  and  April  [14]  have  shown  that  when  an  image 
intensity  is  logarithmically  transformed,  the  speckle  noise 
is  approximately  Gaussian  additive  noise,  and  it  tends  to 
a  normal  probability  much  faster  than  the  intensity  distri¬ 
bution.  Thus,  by  taking  the  logarithm  of  a  speckle  image 
we  get: 

log  /(a:,  y)  =  log  S(x,  y)  +  log  y{x,  y),  ( x,y)eZ 2  (1) 

where  by  I(x,y)  we  denoted  a  noisy  observation  (i.e.  the 
recorded  ultrasound  image)  of  the  two-dimensional  func¬ 
tion  S(x,y)  (i.e.  the  noise- free  image  that  has  to  be  recov¬ 
ered)  and  by  y{x,  y)  the  corrupting  multiplicative  speckle 
noise. 
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Fig.  1.  Examples  of  APD  plots  of  wavelet  subband  coefficients  cor¬ 
responding  to  four  different  ultrasound  images.  In  each  example 
the  empirical  APD  (dotted  line)  is  fitted  with  a  SotS  model  (solid 
line)  and  the  “generalized”  Laplacian  function  (dashed  line). 


The  wavelet  transform  is  a  linear  operation.  Con¬ 
sequently,  after  applying  the  dyadic  wavelet  transform 
(DWT)  to  (1),  we  get  sets  of  noisy  wavelet  coefficients  writ¬ 
ten  as  the  sum  of  the  transformations  of  the  signal  and  of 
the  noise: 

dj,k  =  sj,k  +  €j,ki  (2) 

where  k  =  0,...,2J+-i  —  1  and  —  1  <  j  <  —J  refer  to  the 
decomposition  level  or  scale  and  i  =  1,  2,  3  refers  to  the 
three  spatial  orientations  (see  e.g.  [15]). 

The  signal  components  of  the  wavelet  decomposition 
in  various  scales  are  modeled  as  symmetric  alpha-stable 
(SaS)  processes.  The  SaS  distribution  is  best  defined  by 
its  characteristic  function: 

<j>(u)  =  exp(^w-7|a;|“),  (3) 

where  a  is  the  characteristic  exponent,  taking  values  0  < 
a  <2 ,6  (—oo  <  6  <  oo )  is  the  location  parameter,  and  7 
(7  >  0)  is  the  dispersion  of  the  distribution. 

The  SaS  model  is  suitable  for  describing  signals  that 
have  highly  non-Gaussian  statistics  and  its  parameters  can 
be  estimated  from  noisy  observations.  For  an  experimental 
proof  of  the  appropriateness  of  the  use  of  SaS  models  in 
the  context  of  ultrasound  images  we  refer  the  reader  to  [1] . 
For  the  purpose  of  this  paper  we  only  show  in  Fig.  1  sev¬ 
eral  examples  of  amplitude  probability  density  (APD)  plots 
corresponding  to  subband  representations  of  different  ul¬ 
trasound  images.  Each  data  set  is  modeled  using  both 
the  SaS  family  and  the  “generalized”  Laplacian  density 
function  [7],  [15].  The  plots  prove  that  the  class  of  SaS 
distributions  is  superior  because  it  provides  a  better  fit  to 
both  the  mode  and  the  tails  of  the  empirical  density  of  the 
actual  data. 


III.  Design  of  a  MAP  processor  for  speckle 

MITIGATION 


Motivated  by  the  above  considerations,  we  model  the 
signal  component  of  the  wavelet  coefficients  using  a  two- 
parameter  SaS  distribution1,  while  we  use  a  zero-mean 
Gaussian  model  for  the  noise  component.  Also,  we  consider 
the  signal  and  noise  components  to  be  independent. 

In  a  Bayesian  framework,  referring  to  (2),  dj tk,  and 

£j'k  are  considered  as  samples  of  the  random  variables  d, 
s,  and  £,  respectively.  Our  goal  is  to  find  the  Bayes  risk 
estimator  s  that  minimizes  the  conditional  risk,  which  is 
the  loss  averaged  over  the  conditional  distribution  of  s, 
given  the  noisy  observation,  d: 


s(d)  =  argrnin  J  L[s,s(d)]Ps\d{s\d)  ds 

Selecting  the  uniform  cost  function: 

rrs  g/jM  =  f  °>  for  |s-s|  <  e 
^  (  1,  otherwise 

the  MAP  estimator  can  be  easily  derived  as  being: 
s(d)  =  argmaxPs|d(s  | d) 


(4) 

(5) 

(6) 


It  is  important  to  underline  at  this  point  that  under  the  loss 
function  in  (5),  expression  (4)  is  well  defined  for  all  SaS 
random  variables  (with  characteristic  exponent  a  taking 
values  in  the  whole  range  0  <  a  <  2). 

Bayes’  theorem  gives  the  a  posteriori  PDF  of  s  based  on 
the  measured  data: 


Ps\d(s\d)  = 


Pd\s(d\s)  Ps(s) 

Pd(d ) 


(7) 


where  Ps(s)  is  the  prior  PDF  of  the  signal  component  of 
the  measurements  and  Pd|s(^ls)  is  the  likelihood  function. 
Substituting  (7)  in  (6),  we  get: 
s(d)  =  argmaxP<ns(d  |s)  Ps(s)  =  argma xPj(d  —  s)Ps(s ) 
=  arg  max  Pc  (£)P,(s)  (8) 


Fig.  2  depicts  the  numerically  computed  MAP  input- 
output  curves  for  five  different  values  of  the  signal  char¬ 
acteristic  exponent,  a,  namely,  a  =  2  (Gaussian  data), 
a  =  1.95  (slightly  non-Gaussian  data),  a  =  1.5,  a  =  1, 
and  a  =  0.5  (considerably  heavy-tailed  data).  Apart  from 
the  case  a  =  2,  all  curves  correspond  to  a  nonlinear  “cor¬ 
ing”  operation,  i.e. ,  large-amplitude  observations  are  es¬ 
sentially  preserved  while  small-amplitude  values  are  sup¬ 
pressed.  This  is  expected  since  small  measurement  values 
are  assumed  to  come  from  signal  values  close  to  zero.  On 
inspecting  Fig.  2  it  can  be  observed  that  for  a  given  ratio 
7 /a,  the  amount  of  shrinkage  decreases  as  a  decreases.  The 
intuitive  explanation  for  this  behavior  is  that  the  smaller 
the  value  of  a,  the  heavier  the  tails  of  the  signal  PDF  and 
the  greater  the  probability  that  the  measured  value  is  due 
to  the  signal. 

In  order  to  be  able  to  construct  the  MAP  processor 
in  (8) ,  first  one  should  estimate  the  parameters  of  the  prior 
distributions  of  the  signal  and  noise  components  of  the 
measurements.  Then,  the  parameters  are  used  to  “build” 
the  two  prior  PDFs  P^(£)  and  Ps(s )  and  the  nonlinear  (in 
general)  input-output  relationship  s(d).  To  achieve  this,  we 


1  One  to  the  properties  of  the  DWT,  the  location  parameter  is  in 
this  case  5  =  0. 
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Fig.  2.  MAP  processor  input-output  curves  for  alpha-stable  signal 
(0.5  <  a  <  2)  and  Gaussian  noise  prior  distributions.  The  dash- 
dotted  line  indicates  the  identity  function. 


observe  that  the  PDF  of  the  measured  coefficients  (d)  is  the 
convolution  between  the  PDFs  of  the  signal  (s)  and  noise 
components  (£).  Consequently,  the  associated  characteris¬ 
tic  function  of  the  measurements  is  given  by  the  product 

of  the  characteristic  functions  of  the  signal  and  noise: 

2 

<f>d(w)  =  exp (— 7s|w|as)  •  exp(-y  |w|2)  (9) 

At  this  point,  instead  of  directly  fitting  the  Fourier  trans¬ 
form  of  the  empirical  PDF  of  the  measured  coefficients  with 
the  function  ^(w)  as  we  did  in  [1],  we  observe  that  (9)  im¬ 
plies: 

log  [-(log  |$d(w)|2  +  o"2u;2)]  =  log(27s)  +aslog|w|  (10) 
First,  we  estimate  the  level  of  noise  a  as  in  [6],  then 
we  find  the  parameters  as  and  7S  by  regressing  y  = 
l°g[— (log  |$d(w)|2  +  <t2oj2)]  on  w  =  log  |w|  in  the  model 

yk  =  n  +  a  wk  +  £fc  (11) 

where  y  =  log (27),  ek  denotes  an  error  term,  and  (w*,,  k  = 
1 ,K)  is  an  appropriate  set  of  real  numbers.  We  should 
note  here  that  Koutrouvelis  [10]  used  a  similar  approach  to 
estimate  the  parameters  of  alpha-stable  distributions.  He 
proved  that  its  regression  method  gives  very  good  results 
in  terms  of  consistency,  bias,  and  efficiency. 

IV.  Results 

In  this  section,  we  show  simulation  results  obtained  by 
processing  one  ultrasound  image,  randomly  chosen  from 
our  database.  The  original  image  is  shown  in  Fig.  3(a)  and 
it  represents  an  ultrasound  scan  of  a  fetal  chest. 

In  order  to  obtain  speckle  images,  we  degraded  the  orig¬ 
inal  test  image  by  multiplying  it  with  unit-mean  random 
fields,  as  explained  in  Section  II.  We  controlled  the  cor¬ 
relation  length  of  the  speckle  by  appropriately  setting  the 
size  of  the  kernel  used  to  introduce  correlation  to  the  un¬ 
derlying  Gaussian  noise.  In  practice  uncorrelatedness  of 
the  noise  could  be  achieved  by  decimating  the  image  to 
the  theoretical  resolution  limit  of  the  imaging  device  [13]. 
Thus,  a  short-term  correlation  obtained  with  a  kernel  of 


size  three  was  sufficient  to  model  reality.  We  considered 
three  different  levels  of  simulated  speckle  noise. 

We  compared  the  results  of  our  approach  with  the  clas¬ 
sical  median  filter,  and  wavelet  shrinkage  clenoising  using 
soft  thresholding.  The  soft  thresholding  scheme  was  devel¬ 
oped  using  Daubechies’  Symmlet  8  mother  wavelet  as  sug¬ 
gested  in  [8] ,  while  for  our  algorithm  we  used  the  Symmlet  4 
wavelet.  The  maximum  number  of  wavelet  decompositions 
we  used  was  five.  In  order  to  minimize  the  effect  of  pseudo- 
Gibbs  phenomena,  we  have  embedded  both  wavelet-based 
methods  into  the  cycle  spinning  algorithm  [16]. 

In  order  to  quantify  the  achieved  performance  improve¬ 
ment,  two  measures  were  computed  based  on  the  original 
and  the  denoised  data.  For  quantitative  evaluation,  we 
used  the  signal-to-mean-square-error  (S/MSE)  ratio,  de¬ 
fined  as: 


K  K 

S/MSE  =  101og10(£  S2/J2&  -  ^)2)  (12) 

2=1  2=1 

This  measure  corresponds  to  the  classical  SNR  in  the  case 
of  additive  noise.  Remember  that  in  ultrasound  imaging, 
we  are  interested  in  suppressing  speckle  noise  while  at  the 
same  time  preserving  the  edges  of  the  original  image  that 
often  constitute  features  of  interest  for  diagnosis.  For  this 
reason,  we  also  considered  a  qualitative  measure  for  edge 
preservation  [9]: 


P 


_ T  (AS  -  AS,  AS  -AS) _ 

\Jr(AS-AS,AS-AS)-T(AS-AS,AS-AS) 


_  (13) 

where  AS  and  AS  are  the  high-pass  filtered  versions  of  S 

and  S  respectively,  obtained  with  a  3  x  3-pixel  standard 
approximation  of  the  Laplacian  operator,  and 

K 


T(S1,S2)  =  J2su-S2i.  (14) 


2=1 

The  correlation  measure,  P,  should  be  close  to  unity  for  an 
optimal  effect  of  edge  preservation.  The  results  are  summa¬ 
rized  in  Tables  I  and  II  respectively.  It  can  be  seen  that  our 
proposed  Bayesian  approach  exhibits  the  best  performance 
according  to  both  metrics. 


TABLE  I 

Quantitative  image  enhancement  measures  obtained  using 

THREE  DENOISING  METHODS.  THE  TABULATED  S/MSE  METRIC  IS 
GIVEN  IN  dB. 


Noisy  image 

5.63 

9.67 

16.68 

Median  Filtering 

11.07 

14.24 

17.97 

Soft  Thresholding 

10.95 

14.70 

18.38 

Bayesian  Denoising 

11.66 

15.48 

19.93 

TABLE  II 

Qualitative  image  enhancement  measures  obtained  using 

THREE  DENOISING  METHODS.  VALUES  OF  (3  CLOSE  TO  UNITY  DENOTE 
OPTIMAL  EDGE  PRESERVATION  PERFORMANCE. 


Noisy  image 

0.2577 

0.3930 

0.6933 

Median  Filtering 

0.1989 

0.4230 

0.5310 

Soft  Thresholding 

0.2806 

0.4985 

0.7391 

Bayesian  Denoising 

0.3814 

0.5203 

0.7957 

4 


(d)  (e) 

Fig.  3.  Results  of  various  speckle  suppressing  methods,  (a)  Orig¬ 
inal  image,  (b)  Image  degraded  with  simulated  speckle  noise 
( S/MSE  =  9.67 dB).  (c)  Median  filtering,  (d)  Translation- 
invariant  soft  thresholding,  (e)  Bayesian  denoising. 

In  Fig.  3  we  show  for  visual  comparison  a  representative 
result  from  the  processing  of  our  test  image.  Although  it 
achieves  a  good  speckle  suppression  performance,  the  me¬ 
dian  filter  looses  many  of  the  signal  details  and  the  resulting 
image  is  blurred  (Fig.  3(c)).  On  the  other  hand,  the  image 
processed  by  soft  thresholding  is  oversmoothed  (Fig.  3(d)). 
It  seems  that  the  Bayesian  processor  performs  like  a  feature 
detector,  retaining  the  features  that  are  clearly  distinguish¬ 
able  in  the  speckled  data  but  cutting  out  anything  which 
is  assumed  to  be  constituted  by  noise  (Fig.  3(d)). 

V.  Discussion 

We  introduced  a  novel  technique  for  speckle  noise  re¬ 
moval  using  nonlinear  processing  of  wavelet  coefficients. 
The  proposed  processor  is  based  on  solid  statistical  theory, 
and  it  does  not  depend  on  the  use  of  any  ad  hoc  parameter. 

Our  algorithm  was  tested  and  found  to  be  more  effective 
than  thresholding  methods,  which  do  not  allow  for  an  exact 


matching  of  the  signal  and  noise  distributions  at  differ¬ 
ent  scales  and  orientations.  The  method  proposed  in  Sec¬ 
tion  III  for  choosing  the  “coring”  nonlinearity  could  be  thus 
considered  as  a  principled  way  of  shrinking  noisy  data,  rely¬ 
ing  on  the  true  statistics  of  the  signal  and  noise  wavelet  co¬ 
efficient.  For  example,  the  curve  corresponding  to  a  =  0.5 
in  Fig.  2  mimics  quite  accurately  the  flavor  of  a  hard  thresh¬ 
olding  operator. 

Finally,  we  note  that  our  algorithm  could  be  easily 
adapted  for  the  purpose  of  denoising  other  types  of  biomed¬ 
ical  images  where  the  noise  can  be  (eventually  after  an 
appropriate  transformation)  modeled  as  additive  Gaussian 
and  signal-independent. 
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